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1. Introduction 

The purpose of the Rosetta mission is the in situ analysis of a cometary nucleus using both 
remote sensing equipment and scientific instruments delivered to the comet surface by a lander and 
transmitting measurement data to the comet-orbiting probe. Following a tour of planets including 
one Mars swing-by and three Earth swing-bys, the Rosetta probe is scheduled to rendezvous with 
comet 67P/Churyumov-Gerasimenko in May 2014. The mission poses various flight dynamics 
challenges, both in terms of parameter estimation and maneuver planning. 

Along with spacecraft parameters, the comet’s position, velocity, attitude, angular velocity, inertia 
tensor and gravitational field need to be estimated. The measurements on which the estimation 
process is based are ground-based measurements (range and Doppler) yielding information on the 
heliocentric spacecraft state and images taken by an on-board camera yielding information on the 
comet state relative to the spacecraft. The image-based navigation depends on the identification of 
cometary landmarks (whose body coordinates also need to be estimated in the process). The paper 
will describe the estimation process involved, focusing on the phase when, after orbit insertion, the 
task arises to estimate the cometary rotational motion from camera images on which individual 
landmarks begin to become identifiable. 


2. Scenario Studied 

We choose a space-fixed reference system and denote by £(£) = (£i(£), ^(t), ^(t)) E R 3 the 
coordinate representation of the position vector of the comet’s barycenter at time t with respect 
to this system. Next we choose a coordinate system rigidly attached to the comet, denote by 
gi(t), g 2 (t), gs(t) E R 3 the coordinate representations of these body-fixed directions with respect to 
the reference system and call the matrix g(t) = (gi(t) \ g 2 (t) \ <73 (i)) E SO(3), whose columns are 
the directions gi(t), the comet attitude at time t. It is important to realize that this body-fixed 
system can be arbitrarily chosen. (From a physical point of view it would be most convenient to 
choose gi,g 2 ,g 3 as the directions of the comet’s principal axes, but these principal axes are not 
known initially. From a practical point of view one will rather choose a geometrically defined body- 
system whose direction vectors are defined in terms of landmarks on the comet surface which can be 
identified on camera images.) Finally, we denote by uo(t) E R 3 the body-referenced angular velocity 
vector and by I E R 3x3 the (time-invariant) body-referenced inertia tensor of the spacecraft. To 
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summarize, we introduce the following quantities: 


£(£) = position of comet’s barycenter at time t; 
g(t ) = comet attitude at time t; 

Lo(t) = body-referenced angular velocity vector of comet; 
I = body-referenced inertia tensor of comet. 


( 1 ) 


Moreover, we associate with each vector wel 3 the cross-product operator 
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( 2 ) 


which is characterized by the equation L(oo)v = lo x v for all «Gl 3 . Under the assumption that 
there are no external torques acting on the comet (for example due to outgassing), the rotational 
motion of the comet is governed by the equations 


g = gL(ui) and ICo = ( Ilo ) xw. (3) 

In general, closed-form solutions to Euler’s equation ICo = (Ilo) x lo are not available in terms of 
elementary functions. This is different, however, in the special case that lo points in a space-fixed 
direction which means that Lo(t) = u(t)a where u is a scalar function and where a is a fixed vector. 
In this case the equation ICo = (Ilo) x lo becomes ii(Ia) = u 2 (Ia) x a which is only possible if 
u = 0, i. e., if lo is constant. Then necessarily Ilo x lo = 0 so that lo is an eigenvector of I, i. e., 
is aligned with a principal axis direction. For constant lo the equation g = gL(ui) has the explicit 
solution g(t) = go exp (tL(co)) where go = g(0) is the attitude at the reference time t = 0. 

We now proceed to derive the measurement equations. We assume that a landmark located on the 
comet surface (such as a crater) can be identified on a CCD image taken by the on-board camera. 
If b E R 3 is the body-referenced position of the landmark, then the space-referenced position of 
this landmark at time t is I(t) = £( t ) + bigi(t) + & 2 fl r 2 (^) + = £,(t) + g(t)b. We identify 

the spacecraft position x(t) with the position of the optical center of the on-board camera, denote 
by / the focal width of this camera and by c(t) = (ci(t) | 02 (f) \ 03 (f)) G SO (3) the camera 
attitude (where it is assumed that 03 points in the direction of the optical axis). To summarize, 
we introduce the following quantities: 


b = body-referenced landmark position, 

I(t) = £( t ) + g(t)b = space-referenced landmark position at time t, 
x(t) = spacecraft position (optical center) at time t, (4) 

c(t) = camera attitude at time t, 
f = focal width of camera. 

The image point is given as the point of intersection between the ray from the landmark through 
the optical center and the image plane, which means that there are real numbers A > 0 and «,veR 
(where u and v denote the horizontal and vertical image coordinates) such that I + X(x — I) = 
x + uci + VC2 — fc 3 , as is shown in Figure 1. 
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The equation £ + X(x — £) = x + uc\ + VC2 — fc 3 can be rewritten in the form 

(X-l)(x-£) = uci+vc2-fc 3 . (5) 


Taking the inner product with c 3 we find that (A — l){x — £, c 3 ) = — / and hence that 

A-l = ~ f 

(x - £, c 3 ) ' 

Plugging (6) back into (5) and taking the inner product with c\ and C 2 we find that 


( 6 ) 
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( 7 ) 


For the purposes of our present analysis we may assume that the spacecraft position x is known from 
orbit determination based on ground-based measurements (range and Doppler) and that the camera 
attitude c(t ) is known from spacecraft attitude determination. The task is then to estimate from 
the available measurements (7) (and the underlying dynamical equations) the comet’s position, 
velocity, attitude and angular velocity at some reference time and the body-coordinates of the 
landmarks involved. We first study the phase in which the comet can be identified as a point on 
camera images; in this case the comet itself is treated as a single landmark (with 6 = 0), and the 
task is reduced to estimating the comet position and velocity. 


3. Estimation of Comet Position and Velocity 

We assume that images are taken during a time interval which is short enough to replace the 
spacecraft and the comet orbits by their tangent lines. Assuming thus a linear regime, a moment’s 
thought shows that it is not possible to determine the comet state from a sequence of pictures 
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alone, as the spacecraft-comet vector and the speed of the comet are uniquely determined only up 
to a common scaling factor. 



Fig. 2: Impossibility of determining the comet state simply from camera images. 

What can be done is performing a maneuver (more or less in the blind) which changes the direction 
of the spacecraft motion. In this case the comet state can be determined from four images of which 
at least one each must be taken before and after the maneuver. To derive the necessary formulas, 
let us fix the reference time t = 0 at the time of the maneuver and let us introduce the following 
data: 

xq '■= spacecraft position at maneuver time; 
v : = spacecraft velocity (assumed constant during observation period); 

:= comet position at maneuver time; 

u := comet velocity (assumed constant during observation period);. (8) 

d := £;o — Xq = vector from spacecraft to comet at maneuver time; 
w := u — v = velocity of comet relative to spacecraft before maneuver); 

Av : = velocity change caused by maneuver. 

Then at any time s < 0 the comet and spacecraft positions are £o + su and xo + sv, respectively, 
whereas at any time t > 0 these positions are £o + tu and xq + t(v + Av ), respectively. Thus 
the information to be gleaned from camera images are, for each s < 0, the unit vector e s in the 
direction of (£o + su) — (xq + sv) = d + sw, and, for each t > 0, the unit vector e* in the direction 
of (£ 0 + tu) ~ (xo + tv + tAv) = d + t(w — Av), i.e., the unit vectors 

d + sw d + t(w — Av) 

e s := n and e* := 77- ; — -7 . (9) 

\\d + su> \\d + t(w — Av)\\ 
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Fig. 3: Effect of orbit maneuver. 



Fig. 4: Determination of comet state using an orbit maneuver. 
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Let us assume that images of the comet are obtained at times si < s 2 < S 3 < 0 < t. With the 
notation introduced above, let us write := e Si for the unit vector from the spacecraft to the 
comet at time Sj (where i = 1,2,3). Knowledge of e\ and e 2 fixes the plane in which d and w he. 
Let us write di := \\d + Sjw|| for the spacecraft-comet distance at time Si and let us denote by ipij 
the angle from a to e-q then the Theorem of Sines yields 


sin(<£i 2 ) _ sin(7r - a- <^i 2 ) _ sin(a + tp 12 ) 

(s 2 -si)M| ~~ di di a11 

sin(y?i 3 ) _ sin(-7r -a - ip 13 ) _ sin(o; + ip 13 ) 

(s 3 -si)IMI dx dx 

Taking quotients, we find that 

sin(<£i 2 ) • (s 3 -si) _ sin(a + (^i 2 ) _ tan(ct) cos(<^i 2 ) + sin(<^i 2 ) 
sin(<^i 3 ) • (s 2 -si) sin(a + (^i 3 ) tan(a) cos(<^i 3 ) + sin(<^i 3 ) 


( 10 ) 


( 11 ) 


from which tan(ct) can be determined (and hence also a with a single ambiguity to be resolved). 
Once a is known, (10) relates the speed ||w|| with the distance dx, say ||w;|| = Adi with a known 
factor A. Then the Theorem of Cosines yields 

||d || 2 = d± + s±\\w\\ 2 — 2dx\sx\\\w\\cos(a) 

= d\ + s 2 d 2 A 2 — 2d 2 |si|A cos(a) ( 12 ) 

= d\ (l + sfA 2 — 2|s 1 |Acos(o')) 

so that ||d|| = 0di with a known factor 0. Moreover, for k = 2, 3 the Theorem of Sines yields 


sin(a) sin(7r — a — tpxk) 

dk dx 

dk = Okdx where Ok : = 


sin(o; + y?ifc) 
dx 
sin(a) 

sin(o; + <£i fc ) ' 


and hence 


(13) 


We now make use of the equations d + siw = diei where % = 1,2, 3. We take two of these equations, 
say with indices % and j; subtracting these equations and also subtracting the s^-fold of the first 
equation from the Sj-fold of the second equation, we find that 


w = 


di&i dj&j 0, Cj Oj&j dj Cj Sj d r c r s , 0j Cj SjOjCj 

= di and a = = ai 


Si~Sj 


Si~Sj 


Si~Sj 


Si~Sj 


(14) 


This shows that the measurements taken before the maneuver enable us to determine w and d up 
to a common scaling factor dx- (In particular, the unit vectors = d/||d|| and e w = w/||w;|| can be 
determined from the measurements.) The scaling factor dx can now be determined from one single 
measurement taken after the maneuver; if et is the unit vector from the spacecraft to the comet 
at some time t > 0, we have d + t(w — Av) = fie t for some factor /i. Since d = \\d\\ed = ©die^ and 
w = ||u>||e M = Xdxe w this means 0die^ + t(Xdxe w — Av) = fie t , i.e., 


dxiQea + tXe w ) — fie t = tAv 


(15) 
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where the coefficients d\ and p are unknown whereas everything else is known. Thus d\ and p 
can be found by decomposing tAv into its components relative to the known vectors 0e^ + tXe w 
and et- (The velocity change Av can be assumed to be known from spacecraft orbit determination 
using ground-based measurements after the maneuver.) Hence from now on we can assume that 
the comet’s position and velocity are known and concentrate on estimating the comet’s rotational 
parameters. 


4. Estimation of Cometary Rotation Parameters 

Here we consider the phase in which the spacecraft is in orbit around the comet and the first 
individual landmarks can be identified on CCD images; we try to estimate the rotational motion 
of the comet and the locations of the observed landmarks. We make the simplifying assumption 
that the comet rotates with a constant spin rate co about an axis given by a space-fixed direction d; 
this assumption is legitimate at least for short observation periods. Choosing a body-fixed system 
whose third direction coincides with the (unknown) spin axis direction and denoting by go the 
comet attitude at initial time to = 0, the attitude at any time t is given by 


9(t) 


go 


cos (cot) 
sin(a;£) 


0 


sin(a;£) 0 
cos (cot) 0 
0 1 


(16) 


Hence if a landmark is located at radius r, longitude p und latitude 9 with respect to this body-fixed 
system, then its space-referenced position is given by 



cos 9 cos p 


cos 9 cos(p + cot ) 

m = m+r-g(t) 

cos 9 sin p 

= Z(t)+r-go 

cos 9 sin(<£ + cot) 


sin 9 


sin 9 


(17) 


Let us denote by n the unit vector from the optical center to the landmark, by e the unit vector 
from the comet center to the landmark and by d := £ — x the vector from the spacecraft to the 
comet center; note that both d and n are known, n being given by 


-uci - vc 2 + fc 3 
\/u 2 + V 2 + f 2 


(18) 


Then there is a number A > 0 such that re = — d+Xn. Taking norms on both sides of this equation, 
we find r 2 = ||d|| 2 — 2A(n, d) + A 2 and hence 


Ai ,2 = (n,d) ± y (n,d) 2 - \<l\ - + r 2 . 


(19) 


Since (n,d) > 0 and r < ||d||, both solutions Aj are positive; it is geometrically clear that, if the 
comet shape is not too different from that of a sphere, the smaller of the two values A j is the sought 
solution. Hence A = (n, d) — ( n , d) 2 — ||d|| 2 + r 2 so that the equation re = — d + An becomes 

e = - d + (n, d)n — \j (n, d) 2 — ||d|| 2 -|-r 2 • nj . (20) 


7 




Fig. 5: Determination of the unit vector from the comet center to the landmark 

Let us now assume that the comet can be approximated by a sphere of radius R where R is at 
least approximately known. Then r = R. and the upshot of (20) is that the unit vector e(t) can be 
computed from a measurement taken at time t. This observation will now be exploited to determine 
the spin rate cu, the comet attitude go at the reference time and the landmark coordinates 9 and 
ip from given landmark vectors ej = e(ti) resulting from measurements at times t 4 . 

Method 1: One landmark from three images. Take three measurements on which the 
same landmark can be identified and let ej, ej and e k be the associated landmark unit vectors. 
The comet’s rotation axis d must be perpendicular to both ej — ej and e k — c.j . hence is given by 

j = (ej-ej) x (e k - ej ) 

\\ {ej - e 4 ) x (e k - ej )\\ ■ [ ^ 



The landmark latitude 9 must then satisfy sin 2 9 = cos 2 (90° — 9) = (d, e^) 2 so that 

sin 9 = ±(d,ei ) (22) 


ej - (ej, d)d 
ei - {ej, d)d\\ ’ 

8 


for all i. Moreover, writing 


( 23 ) 


and denoting by ipij the angle by which the landmark was rotated about the axis of rotation, we 
have 

(e!,e}) = cos < fij = cos (u(tj - U)) (24) 

which determines the spin rate to. Consequently, the comet attitude at time ti is given by the 
column vectors 

d & ' 

gs = d, g 2 = -rrz ^ and gi = g 2 x g 3 . (25) 

It is clear that the landmark longitude cannot be estimated from the available measurements; 
however, we can assume that the first axis in the body-fixed system is pointing in the direction of 
the landmark longitude so that ip = 0 (i.e. , we define zero longitude by the landmark location). 

Method 2: General equations. Let us assume that a number of landmarks can be identified 
on a sequence of images. Let be the unit vector from the comet center to the i-th landmark 
at the time the k - th image is taken. Dividing (17) by r = R we find that 


e b,fc) 


Oo 


COS Oi COS ( pi + LOtk ) 

cos Oi sin(<£j + Lotk) 
sin Oi 


(26) 


where tk is the time (counted from the epoch) at which the k- th picture is taken. Since go is a 
rotation matrix and hence preserves inner products, we find from (26) that 

= cos Oi cos Oj cos(pi — tpj + Lo(tk — ti)) + sin Oi sin Oj (27) 

so that the comet attitude go is decoupled from the other unknowns. (We note that it is not possible 
to determine all longitudes pi from these equations, since only the differences pi — <pj occur; this 
shows again that, within our present considerations, zero longitude can be defined arbitrarily for 
one of the landmarks, and one can only hope to determine the longitude separations from all other 
landmarks to the chosen one.) Let us see which conclusions can be drawn from (27) in special 
cases. 

Method 2a: One single landmark from several images. Considering only one single 
landmark (say the i-th), we find that 


(e ( -*’ fc \ e^ 1 ’^) = cos 2 Oi cos (cotk — coti) + sin 2 Oi 
= cos (cotk — coti) + (l — cos (cotk — uti)) sin 2 Oi 


Eliminating Oi via 


sin 2 Oi 


( e b,fc), gbh)) _ cos (cotk ~ uti) 
1 — cos (cotk — coti) 


we obtain the compatibility conditions 


( e (*,fc), e ( l ,£)) _ cos (cotk — tut i) 
1 — cos (cotk — coti) 


(g(*’ K ), g(*’ L )) — COS (cotK — OOtL,) 
1 — COS (cotK — ooti) 


(28) 


(29) 


(30) 


which must hold for all pairs (k,£) and ( K,L ). Note that these can be rewritten in the form 
akt cos(loA kl ) — a K L cos (lo A k i ) +a K L - a k i = 0 where A k t '■= tk~U and a k t := 1 - e ( *’ £) ); 

hence lo can be found as a root of functions of the form F(lo) := a cos(Ao;) — b cos {Boo) +b — a with 
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known constants a, b, A, B. Thus the following algorithm can be used to reconstruct the unknown 
parameters from one single landmark which is seen on at least three different images: We first 
determine the spin rate lo from (30) (possibly with ambiguities), then the landmark latitude 0i 
from (29) (with a single ambiguity) and finally the comet attitude go from (26) (to be solved in 
the least-squares sense after letting ^ := 0 ). 

Method 2b: Several landmarks on one single image. Let us see which information can 
be inferred from one single image (say the k-th). Letting £ = k in (27), we find that 

(e^ ,k \e^ l,k ^) = cos 0i cos 0j cos((pi — (fij) + smOismOj. (31) 


If the landmark latitudes 0i are already known (for example from the method described before), 
we can simply use (31) to determine the longitude separations via 


COS (<Pi - (fij) 


(e OT) , e 0 , fc ) ) — sin 0i sin 0j 
COS 0i COS 0j 


(32) 


If (32) is used without a priori information then, if s landmarks can be identified, this constitutes 
a system of s(s — l )/2 equations for the 2 s — 1 unknowns 0±, . . . , 0 S and <^ 2 — <Pi, ■ ■ ■ , ‘fs ~ <Pi, which 
can be solved if s(s — l)/2 > 2s — 1 (which is the case for s > 5). 


Method 2c: Two landmarks from two images. Let us consider the case that two different 
landmarks can be identified on two different images. Writing At := t 2 — and A ip : = tp 2 — ‘Pi, 
the nontrivial equations of the form (27) (i.e., the ones for which (i, k ) 7 ^ (j,£)) take the following 
forms: 



e* 1 ' 2 *) 


e* 2 ' 1 )) 


e< 2 ' 2 >) 


e* 2 ' 1 )) 

e* 1 ' 2 ), 

e< 2 ’ 2 )) 

e* 2 ' 1 ), 

e< 2 ’ 2 )) 


cos 2 0i cos(coAt) + sin 2 0i; 

cos 0i cos 02 cos(A^) + sin 0i sin 0 2 ; 

cos 0i cos 02 cos(A<^ + a; At) + sin 0i sin 0 2 ; 

cos 0i cos 02 cos(A<^ — a; At) + sin 0i sin 0 2 ; 

cos 0i cos 02 cos(A^) + sin 0i sin 0 2 ; 

cos 2 02 cos(a;At) + sin 2 0 2 . 


(33) 


A comparison of the second and the fifth equation yields the compatibility condition 


(ei 1,1 ) , e( 2,1 )) 


(et 1 ’ 2 ),^ 2 ’ 2 )); 


(34) 


since (33) consists of six equations for four unknowns there must be (at least) one other compati- 
bility condition. From the first and the last equation we find that 


cos 2 0i 


l-<e(i-i),e(i- 2 )> 

— ; 7 — — an d 

1 — cos(LoAt) 


cos 2 02 


l_( e (2,l), e (2,2)^ 

1 — cos(LoAt) ’ 


(35) 


since the landmark latitudes must satisfy — tt/ 2 < 0i < tt/2 and hence have a positive cosine, roots 
can be extracted unambiguously in (35) so that 


cos 0i 


l-(e( 1 . 1 ), e ( 1 .2)) /l-(e( 2 d), e ( 2 - 2 )) 

— — t — - — and cos u 2 — \ — 7 — 7 — r — 

1 — cos(o;A t) y 1 — cos(o;A t) 


(36) 
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Plugging in (33) the second into the third and the hfth into the fourth equation, we find that 


(e*- 1,1 ^, e ( ' 2,2 ' ) ) = cos cos 02 cos(A<^ + coAt) + (e*- 1,1 ^, e*' 2,1 -*) — cos0i cos 02 cos(A^), 
(e*- 1,2 ^, e ( - 2 ’ 1 ' ) ) = cos 0i cos 02 cos(A<^ — coAt) + {e^ 1,2 \ e*' 2 ’ 2 -*) — cos0i cos 02 cos(A^) 


and hence that 


=(1,1) A 2 , 2 ) _ „(2,1)\ _ 


gi 2 , 1 )) _ cos 0i cos 02 • (cos(A<^ + coAt) — COs(A(^)) , 


(e*- 1,2 ^, e 1 - 2 ’ 1 ) — e ( ' 2,2 ' ) ) = cos 0i cos 02 • (cos(A<^ — coAt) — cos(A^)) . 

Plugging (36) into (38), we obtain 

(gib 1 ), e ( 2 ’ 2 ) — e^ 2,1 )) cos(A<^ + coAt) — cos(A< p) 

y/1 - (ei 1 ’ 1 ), ePAiy v / T^~( 6 ( ' 2 ’ 1 ), e( 2 ’ 2 )) ~~ 1 - cos(coAf) 

(g( 1,2 ), e^ 2,1 ) — e^ 2,2 )) cos(A<^ — coAt) — cos(A^) 

y/1 - (ei 1 ’ 1 ), ePAiy v / T^( 6 ( ' 2 ’ 1 ), e( 2 ’ 2 )) ~~ 1 - cos(o>Af) 

Addition of the last two equations yields U + V = — 2cos(A<^), i.e., 

1 (g(l’ 2 ) — g(l’l), g( 2 ’ 2 ) — g(2,l)\ 

COS(A^) = - 


U : = 
V : = 


( 37 ) 


(38) 


(39) 


2 _ (ei 1 ’ 1 ), ei 1 ’ 2 )) y/l - <e( 2 A), e ( 2 , 2 )) ' 

Once A ip is known, we can determine coAt (and hence co) from any of the two equations in 
which can be rewritten in the form 

cos(Aip) cos(coAt) — sm(Aip) sm(coAt) — cos(Aip) = U — U cos(coAt), 
cos(Aip) cos(coAt) + sm(Aip) sm(coAt) — cos(Aip) = V — V cos(coAt). 


(40) 
(39), 

(41) 


A particularly elegant solution (which yields coAt directly in terms of the measurement data, 
without a need to determine A p beforehand) is obtained by rewriting (41) in the form 


cos(coAt) — 1 — sin(a;At) 


cos(A^) 

cos(coAt) — 1 sin(a;At) 


sin(A<^) 


which, upon inversion, results in 


(l — cos(coAt)) 


U 

V 


cos(A^) 

-1 

sin(a;At) sin(a;At) 

~U~ 

sin(A<^) 

2 sin(a; At) 

1 — cos(coAt) cos(coAt) — 1 

V 


i.e., 


2 sin(a;At) 


cos(A^) 


( U + V ) sin(a;At) 

sin(A<^) 


( U — V) (l — cos(a;At)) 


( 42 ) 


( 43 ) 


( 44 ) 


Taking norms on both sides yields 4sin 2 (a;At) = (U + V ) 2 sin 2 (a;At) + (U — V) 2 (l — cos(coAt )) 2 
which can be rewritten as (4 — (U + V ) 2 ) (l — cos 2 (a;At)) = (U — V) 2 (l — cos(coAt)) 2 . Dividing 
this equation by 1 — cos(a;At) yields (4 — (U + F) 2 ) (l + cos(a;At)) = (U — 14) 2 (l — cos(a;At)), 
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which results in (4 + (U — V ) 2 — (U + V ) 2 ) cos(o;At) = (U — V ) 2 + (U + V ) 2 — 4. This can be 
rewritten in the form 


cos(LoAt) 


2 (U 2 + V 2 ) -4 

4-4UV" 


\{U 2 + V 2 ) - 1 

1 - 171/ 


(45) 


Subsequently, the latitudes 9 1 and 02 can be obtained from (36). 

Method 2d: Determination of longitude separations. If the landmark latitudes 9i and 
9j and the spin rate to are already known, we can rewrite (27) in the form 


COS (ifii - (fij + Lo(t k - tl)) 


(£^’ fc ),£^"’ £ )) — sin0isin0j 
cos 9i cos 9j 


(46) 


to determine the longitude separation <pi — ipj if landmark i can be seen on the k - th image and 
landmark j on the Gth image. 


5. Averaging Initial Estimates 

In the previous paragraph it was described how initial estimates for the estimation parameters can 
be obtained under simplifying assumptions. However, different methods and different images used 
may lead to different raw estimates, and to obtain an initial estimate which is good enough to ensure 
convergence of the subsequent estimation process, it is necessary to form an “average” of the raw 
estimates obtained. This leads to a general problem: Given elements p±, . . . ,p n on a manifold M, 
what is the “average value” of these elements? (In our case the manifold M is the three-dimensional 
rotation group, and the elements to be averaged are rotation matrices representing different raw 
estimates of the comet attitude. Similarly, in spin axis attitude determination, the manifold M 
is the unit sphere in three-dimensional space, and the elements to be averaged are unit vectors 
representing different raw estimates of the spin axis direction of a spacecraft.) Average values 
of elements pi, . . . ,p n on a Riemannian manifold were first discussed by Maurice Frechet* who 
defined them as those elements £ which minimize d(^,Pi) 2 where d(p,q) is the Riemannian 

distance between points p and q, i.e., the length of the shortest curve joining p and q. In most 
practical cases such a mean can be found by an iterative scheme whereby an initial guess £ for the 
average of p±, ... ,p n leads to an improved guess £ via the general formula 


f := exp ? ^XJlog ? (pi)j (47) 

where exp^ denotes the exponential function at a point £ (defined by exp^(u) = a( 1) where a is 
the unique geodesic with a(0) = £ and d(0) = v) and where log^ is its inverse. 


* L’integrale abstraite d’une fonction abstraite d’une variable abstraite et son application a la 
moyenne d’un element aleatoire de nature quelconque, Revue Scientifique 1944, pp. 483-512. - Les 
elements aleatoires de nature quelconques dans un espace distancie, Annales de I’Institut Henri 
Poincare 10 (1948), pp. 215-310. 
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Fig. 7: Improving an approximate mean £ of points p ]:, .... p n on a Riemannian manifold M. 

Examples, (a) If M = S n_1 := { x G R n | ||x|| = 1} is the unit sphere in ^.-dimensional 
space, then the exponential function and the logarithm are given by 


exp,(r) = cos(||w||)C + sin(||w| 


. . , arccos ((»,£)) i , .. . 

log € (p) = 


(48) 


V 1 - (P’0 2 


(b) If M = G is a Lie group (such as the rotation group), then everything can be reduced to 
a neighborhood of the identity element e; i.e. , the iteration scheme takes the form 


1 


? = ? ex Pt ( - 1r * 


(49) 


i= 1 


(c) If M = SO(3) := {g G R 3x3 | g T g = 1} is the rotation group in three-dimensional space 
(with the identity matrix 1 playing the role of the neutral element e), then the Rodrigues formula 
yields 


( T , \\ /n in. 1 — cos ( \\lo ) sin( , 

exp e (L(o;)) = cos(||o;||)H frTTT ^ cu (S H n _ n L(uj), 


a; 


\0J\ 


arccos((tr(#) - l)/2) T 

log„o = . =(q — g ). 

\/3 - tr(#) 2 + 2tr(g) 


(50) 
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6. Iterative Improvement of Estimates 

In the previous paragraphs we showed how rough estimates of the comet attitude and angular 
velocity can be obtained. These estimates need to be improved in an iterative scheme whereby 
in each step the measurement equations are linearized about the current estimates and then a 
least-squares minimization of the resulting residuals is performed to update these estimates. For 
the purpose of this paper we assume that both the spacecraft position and the camera attitude are 
known. Then the only uncertain term in the measurement equations 


u 


( L,c i 

<w 


/ ' if'’ Cl | (where L := £ — x) 
\L, c 3 ) 


(51) 


is the landmark location £. A change 8£ = 8L in the nominal value £ results in changes 8u and 8v 
in the expected measurements which are given by 


{6L,c 1 ){L,c 3 )-{L,c 1 ){6L,c 3 ) 

1 <L,C 3 > 2 

{8L,c 2 ){L,c 3 )-{L,c 2 ){8L,c 3 ) 
1 {L, c 3 ) 2 


f 

f 


(C1XC3, ( 8L)xL ) 
(L, c 3 ) 2 

(c 2 xc 3 , ( 8L)xL ) 
(L, c 3 ) 2 


f 

f 


(. Lxc 2 ,8L ) 
(L,c 3 y 
(. Lxc u 8L ) 
(L,c 3 ) 2 


(52) 


Now a change m L = £ — x = £ + gb — x results from changes in g and b. Let uo, go and 
loo be the position, velocity, attitude and angular velocity at a reference time to = 0. While in 
operational software we have to derive ^ and g at the measurement time from these initial data 
by numerically integrating the dynamical equations, we assume here for simplicity’s sake that the 
comet velocity and angular velocity can be treated as constants (which is not too unrealistic for a 
sufficiently short observation interval). Then £(t) = £,o + tu and g(t) = go exp(tL(o;)) and hence 


8£ = 8t ; o + t ■ 8u, 

Sg = (Sgo)exp(tL(Lu)) + t ■ g 0 exp(tL(tj))L(5tj). 


(53) 


Now due to the structure of the rotation group SO (3) as a nonlinear manifold, a change (i. e., 
infinitesimal increment) 8 go in the estimate go takes the form 8 go = goL(A) with a vector A = 
(Ai, A 2 , A 3 ) t G M 3 . Then, writing g = g 0 exp(tL(o;)) , the second equation in (53) becomes 


$9 = 9 oL(A) exp(tL(o;)) + t ■ go exp(tL(o;))L(&<;) 

= SoifAjgjT'g + t ■ gL(Soj) = L(g 0 A)g + t ■ gL(Soj) . 


(54) 


Thus changes in £o, u, go, w and the body-referenced landmark position b result in a change in 
L = £ + gb — x which is given by 


8L = (8£) + (8g)b + g(8b) = (8£ 0 ) +t ■ (8u) + L(g 0 A)gb + t ■ gL(8u)b + g(8b) . (55) 

Now if a sequence of images is taken on which a number of s landmarks (out of which r are 
different) can be identified, then plugging (55) into (52) results in an equation of the form 
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( 56 ) 


Su i ' 


I 

O 

00 ^ 

1 

5-ui 


A 

8u s 

= A 

Su 

8b^ 

- SV s . 


8b W 


where A is a matrix of size (2s) x (12 + 3 r) involving the comet states at the measurement times. 
As it is desired to match, as well as possible, the theoretically expected with the actually obtained 
measurements, this leads to solving the overdetermined linear system 


5£o ' 


Su 


r„, obtained expected " 

Uj - ^ 

A 


...obtained expected 

v 1 — iq 

8u 

= 


or 

1- 


^.obtained ^.expected 

tt ^ tv ^ 



obtained expected 

8b W 


- V s v s 


(57) 


in the least-squares sense. (Here the expected measurement values are obtained by evaluating 
(7) using the current estimates for the estimation parameters.) Once this is done, the computed 
parameter increments 8£ o, Su, A, 8u and 8b^> are used to update the current estimates. In the 
case of the parameter go, this update state is nonlinear, being given by the formula 

9o = £oexp(L(A)). (58) 


go L(A) 



Fig. 8: Nonlinear update step for the estimate of the initial comet attitude. 
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